Nonlinear  Krylov-Secant  Solvers* 

Hector  Kite,  and  Mary  F.  Wheeler 

The  Center  for  Subsurface  Modeling  (CSM) 

The  Institute  for  Computational  Engineering  and  Sciences  (ICES) 
The  University  of  Texas  at  Austin,  Austin,  TX,  78712 


Contents 

1  Introduction  2 

2  Newt  on- Krylov  Framework  3 

3  The  Family  of  Secant  Solvers 

3.1  Broyden’s  Method . 

3.2  The  Nonlinear  Eirola-Nevanlinna  (NEN)  Method 

4  Krylov-Secant  Updates  9 

4.1  Secant  Updates  constrained  to  the  Krylov  subspace .  9 

4.2  Krylov  Secant  vs.  Standard  Secant  Update  .  12 

4.3  The  Nonlinear  Krylov-Eirola-Nevanlinna  Algorithm .  14 

4.4  A  high-order  Newton— Krylov  algorithm  .  15 

5  Implementation  issues  15 

5.1  Preconditioning .  15 

5.2  Globalization  Strategy  and  Forcing  Terms .  16 

6  Numerical  Experiments  17 

7  Conclusions  and  Further  Work  19 

Abstract 

This  report  describes  a  new  family  of  Newton-Krylov  methods  for 
solving  nonlinear  systems  of  equations  arising  from  the  solution  of 
Richards’  equation  and  in  fully  implicit  formulations  in  air-water  sys¬ 
tems.  The  basic  approach  is  to  perform  secant  (Broyden)  updates 
restricted  to  the  Krylov  subspace  generated  by  the  GMRES  iterative 
solver.  This  approach  is  introduced  as  Krylov-secant  methods.  One 
of  the  most  attractive  features  of  these  methods  is  their  performance 
of  sequence  of  rank-one  updates  without  explicitly  recalling  the  com¬ 
putation  or  action  of  the  Jacobian  matrix.  Implications  of  these  up¬ 
dates  in  line-search  globalization  strategies,  computation  of  dynamic 
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tolerances  (forcing  terms)  and  the  use  of  preconditioning  strategies 
are  presented.  Numerical  results  show  improvements  over  traditional 
implement  ations . 


1  Introduction 

It  is  well  known  by  environmental  engineers  that  solving  Richards’s  equa¬ 
tion  can  be  a  challenging  problem  due  to  the  high  nonlinearities  induced  by 
pressure  and  saturation  dependent  coefficients  in  the  temporal  and  spatial 
terms  [21].  Moreover,  when  sufficient  computer  power  is  available  for  solving 
more  complex  air- water  groundwater  applications,  fully  implicit  formulations 
are  preferred  over  IMPES  formulations  since  they  are  unconditionally  stable 
and  allow  the  performance  of  larger  simulation  timesteps  [3,  11,  14].  These 
timesteps  are  numerically  constrained  to  Newton’s  radius  of  local  conver¬ 
gence,  and  the  cost  of  performing  each  one  of  them  may  be  computational 
high. 

On  the  other  hand,  recent  advances  in  Newton-Krylov  methods  [17,  19] 
have  brought  up  new  opportunities  for  environmental  agencies  to  speedup 
the  computation  in  fully  implicit  implementations  [4,  16,  17].  Nowadays,  the 
use  of  Krylov  iterative  methods  [27,  31]  such  as  Orthomin,  QMR,  BiCGSTAB 
and  GMRES  provides  powerful  means  to  implement  inexact  versions  of  New¬ 
ton’s  method  for  solving  large  scale  problems  in  an  efficient  way. 

Despite  these  rapid  advances,  linear  preconditioning  is  still  a  major  issue 
for  fully  implicit  implementations  in  multiphase  flow  due  to  the  coupling  of 
different  physical  variables  [4,  15,  16].  Nevertheless,  it  is  clear  that  improving 
the  nonlinear  convergence  may  result  in  significant  savings  in  CPU  time 
since  each  nonlinear  iteration  typically  involves  dozens  to  hundreds  of  linear 
iterations. 

The  present  work  proposes  an  efficient  way  to  exploit  the  underlying 
Krylov  subspace  information  from  the  linear  solver  in  order  to  accelerate  the 
nonlinear  convergence  of  Newton-Krylov  methods.  The  basic  approach  is 
to  perform  secant  (Broyden)  updates  restricted  to  the  Krylov  subspace  gen¬ 
erated  by  the  GMRES  iterative  solver,  namely,  updates  to  the  Hessenberg 
matrix  resulting  from  the  Arnoldi  factorization.  This  approach  is  introduced 
as  a  new  family  of  Krylov-secant  methods.  One  of  their  most  attractive 
features,  beside  achieving  higher  rates  of  nonlinear  convergence,  is  the  per¬ 
formance  of  sequences  of  rank-one  updates  without  explicitly  recalling  the 
computation  or  action  of  the  Jacobian  matrix.  The  procedure  may  be  de¬ 
scribed  as  the  composition  of  several  chord  steps  for  every  single  Newton 
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step  in  a  predictor-corrector  fashion. 

Therefore,  in  contrast  to  other  fully  implicit  implementations  for  ground- 
water  flow  implementations  this  work  focuses  on  reducing  the  number  of  non¬ 
linear  iterations  at  a  low  computational  cost.  In  fact,  Krylov-secant  methods 
may  be  implemented  as  a  practical  extension  to  the  Newton-Krylov  solver 
framework.  We  discuss  implications  of  these  updates  in  the  formulation  of 
line-search  globalization  strategies,  computation  of  dynamic  tolerances  (forc¬ 
ing  terms)  and  the  use  of  preconditioning  strategies. 

Numerical  results  show  improvements  over  traditional  implementations, 
especially  in  those  cases  where  nonlinearities  increase  due  to  stringent  pres¬ 
sure  and  saturation  changes  throughout  the  simulation. 

Further  potentials  of  the  proposed  method  are  devised  with  the  current 
advances  in  augmenting,  reorthogonalizing  and  shifting  the  Krylov  basis  in¬ 
formation  [5,  23,  26]. 

2  Newton-Krylov  Framework 

We  are  interested  in  solving  the  following  nonlinear  system  of  equations 

*»=0,  (1) 

where  F  :  fl  C  Mn  — >  Rn.  A  variety  of  different  discretization  techniques 
applied  to  nonlinear  partial  differential  equations  give  rise  to  an  algebraic 
system  such  as  (1),  where  u  is  the  unknown  vector  representing  the  nodal 
values  of  the  state  variable. 

In  our  particular  case,  this  nonlinear  system  of  equations  arises  from  dis¬ 
cretizing  either  Richards’  equation  or  the  coupled  partial  differential  equa¬ 
tions  describing  the  air-water  system.  The  nonlinearities  in  both  equations 
are  present  in  the  temporal  term  and  the  diffusive  term,  and,  when  fully 
implicit  schemes  are  employed,  the  transient  and  steady-state  cases  lead  to  a 
nonlinear  algebraic  system  such  as  (1).  Newton’s  method  implies  the  solution 
of  the  Jacobian  system 

J  (it)  5u  =  —  F  (it)  ,  (2) 

where  F  (u)  and  J  (it)  denote  the  evaluation  of  the  function  at  u  £  Rn  and 
its  derivative  at  any  Newton  step,  respectively. 

Due  to  the  high  cost  or  absence  of  derivatives  for  the  explicit  construction 
of  the  Jacobian  J,  the  action  of  this  operator  onto  a  vector  (in  the  form  of 
matrix-vector  products)  may  be  performed  by  finite  differences.  This  is,  in 
fact,  one  of  the  most  appealing  features  of  Newton-Krylov  methods.  These 
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methods  are  also  amenable  to  globalization  by  line-search  or  trust  region 
strategies  and  for  dynamic  control  of  linear  tolerances  (i.e. ,  forcing  terms). 
The  reader  is  referred  to  [19]  to  see  a  state-of-the-art  discussion  of  the  subject. 
Specific  theory  supporting  the  description  of  the  line-search  backtracking 
method  and  the  selection  of  the  forcing  terms  to  control  the  convergence 
of  the  linear  solver  at  each  nonlinear  step  may  be  found  in  [10,  17,  25].  A 
globalized  Newton-Krylov  method  with  a  line-search  backtracking  procedure 
may  be  stated  as  follows: 

Algorithm  2.1  (Newton-Krylov  algorithm  with  forcing  terms  and  line-search 
backtracking) 

1.  Let  u ^  be  an  initial  guess. 

2.  For  k  =  0, 1,  2, . . .  until  convergence  do 

2.1  Choose  e  [0, 1) ,  t  €  (0, 1). 

2.2  Using  some  Krylov  subspace  iterative  method,  compute  a  vector 
s W  satisfying 

J(fcVfc)  =  -F(k)  +r(fc),  (3) 

1 1  7»(k)  II  /  J  \ 

with  1, 11.  ,,  Ci  <  m  b 

2.3  While  +  s(fc))||2  >  [l  -  t  (l  -  r/^)]  do 

2.3.1  Choose  A  to  minimize  a  quadratic  polynomial  over  A,  A  C 
(0, 1)  that  interpolates 

F  (A)  =  ||F(u(fc)  +  As(fc))|2. 

2.3.2  Update  =  A s^. 

2.4  Update  =  1  —  A  ^1  — 

3.  Update  solution  u (fc+1l  =  +  s^k\ 

We  can  make  the  following  observations: 

•  The  loop  at  step  2.3  represents  a  global  backtracking  based  on  the 
Goldstein- Armijo  condition  [7].  In  practice,  the  parameter  t  is  se¬ 
lected  on  the  order  of  10-4.  The  interval  [a,  a]  =  [.1,  .5]  is  typically 
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chosen,  and  the  quadratic  interpolant  p  is  constructed  in  such  a  way 
that  p(0)  =  F( 0)  =  F (fc)  2  ,  p(l)  =  F(l)  =  2and 

p'  (0)  =  F'  (0)  =  2  (f( k\ 

•  The  forcing  term  rfk'1  is  selected  according  to  the  following  criterion: 


where 


r]{k)  =  min  j  r?max,  max  j  r]{k\  [rfk  1}) 


?7(fc)  = 


p(k)  _  -g  !)s(fc— !) 

'  m=m  1 


(4) 

(5) 


The  choice  of  rj ^  reflects  the  agreement  between  F  and  its  linear  model 
at  the  previous  iteration,  and  it  depends  on  how  A  is  chosen.  Thus,  the  linear 
solver  tolerance  is  larger  when  the  Newton  step  is  less  likely  to  be  productive 
and  smaller  when  the  step  is  more  likely  to  lead  to  a  good  approximation. 
This  operation  is  performed  once  the  backtracking  procedure  has  returned 
the  adequate  nonlinear  step  size  to  update  the  solution.  Fig.  1  illustrates 
the  savings  in  computation  using  the  forcing  terms  for  a  fully-implicit  for¬ 
mulation  for  a  3-D  two-phase  problem.  The  staircase  shape  displayed  by  an 
inexact  Newton  method  with  a  fixed  forcing  term  suggests  the  amount  of 
over-computation  in  generating  decreasing  Newton  directions.  In  contrast, 
the  criterion  given  by  (4)  and  (5)  avoids  flat  portions  of  the  curve  resulting 
then  in  a  significant  overall  saving  of  GMRES  iterations  (about  400  itera¬ 
tions)  . 


3  The  Family  of  Secant  Solvers 

Rank-one  updates  for  solving  nonlinear  equation  are  sometimes  referred  to 
as  secant  or  quasi-Newton  methods  since  no  “true”  Newton  equation  is  ever 
formed  throughout  all  iterations.  These  updates  obey  a  secant  condition 
that  must  be  satisfied  by  the  new  Jacobian  approximation.  The  best  of  these 
methods  still  seems  to  be  the  first  one  originally  introduced  by  Broyden  (see 
e-g.,  [7]). 
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Figure  1:  The  use  of  the  forcing  term  criterion  for  dynamic  control  of  linear 
tolerances.  The  solid  line  represents  a  standard  inexact  Newton  implemen¬ 
tation  with  fixed  linear  tolerances  0.1.  The  dotted  line  is  the  inexact  Newton 
implementation  with  the  forcing  term  criterion.  Each  symbol  *  indicates  the 
start  of  a  new  Newton  iteration. 
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3.1  Broyden’s  Method 

Given  u  ~  u*  and  M  ~  J  ( u ) ,  we  can  find  an  approximate  new  Newton  step, 
u+,  by 

u+  =  u  —  A~lF  ( u ) .  (6) 

Broyden’s  method  computes  a  new  A+  by  means  of  the  following  rank- 
one  update 

A+  =A+  1F+  (“)  -  F  (“>  -  Ag]  - ,  (T) 

gls 

that,  whenever  the  approximated  Jacobian  equation  is  solved  exactly,  i.e. 
As  =  —F(u),  reduces  to 


A+ 


A  + 


F+  (uW 
g*s 


for  g^s  ^  0.  The  vector  g  may  be  chosen  in  several  ways.  For  instance,  when 
g  =  s,  we  obtain  the  “good  Broyden’s  update”  and  when  g  =  A1  [i?+  (it)  —  F  (it)], 
we  have  the  “bad  Broyden’s  update”. 

One  of  the  key  features  of  Broyden’s  method  is  to  provide  a  sparse  ap¬ 
proximation  to  the  true  Jacobian  that  may  be  exploited  in  structure.  In  fact, 
there  are  efficient  procedures  to  factorize  such  an  operator  under  low-rank 
modifications  [7]. 

In  the  setting  of  groundwater  flow  problems,  Broyden’s  method  may  be  a 
convenient  option  to  avoid  the  computation  of  several  hundred  thousands  to 
millions  of  entries  in  the  Jacobian  matrix.  However,  caution  must  be  taken 
when  the  updated  operator  A+  becomes  numerically  singular.  In  that  case, 
the  Jacobian  must  be  reconstructed  at  the  new  point  to  resume  the  nonlinear 
iterations. 

This  algorithm  converges  g-super linearly  to  F  ( u *)  if  the  rate  at  which 
the  Jacobian  is  approximated  is  faster  than  the  reduction  in  the  nonlinear 
step  (i.e.,  the  Dennis- More  characterization  [7]). 


3.2  The  Nonlinear  Eirola-Nevanlinna  (NEN)  Method 

The  nonlinear  Eirola-Nevanlinna  (NEN)  was  proposed  by  [34]  as  the  nonlin¬ 
ear  counterpart  to  the  linear  EN  algorithm  [9].  This  method  is  a  general¬ 
ization  of  Broyden’s  method  for  solving  linear/nonlinear  equations  and  has 
received  particular  attention  in  recent  years  [13,  32],  Interestingly  enough, 
the  algorithm  has  connections  with  the  GMRES  algorithm  [31].  The  NEN 
algorithm  may  be  regarded  as  the  composition  of  two  Broyden’s  iterations 
as  the  following  presentation  suggests: 
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Algorithm  3.1  (Nonlinear  EN) 


1.  Give  an  initial  guess  u ^  and  Jacobian  approximation 

2.  For  k  =  0, 1, . . .  until  convergence  do 


2.1 

2.2 

2.3 

2.4 

2.5 


Solve  = 

q(k)  _  ^(fc+i)  _ 

J^(k+l)  _  J^{k)  _j_ 

Solve 


:  — 

F(k). 

(,(fc)_AW,W)(aW)t 

(sC^) ^  ^s(^) 


Update  solution 


Note  that  each  of  the  iterations  involves  the  solution  of  two  linear  sys¬ 
tems.  However,  we  still  maintain  one  rank-one  update  and  one  function 
evaluation  per  step  as  in  Broyden’s  method.  It  can  be  shown  that  the  up¬ 
date  at  step  2.5  can  be  expressed  as 


where 


and 


u(k+i)  =  u(k)  +  jffc) 

=  u(fc)  +  +  0(fc>5<fc\ 


(*)  -  _ 


sv  '  = 


(A(fc))  1  F{k\ 


a(fc)  =  - 


F  (u(fc)  +  s(fc))  , 

g(k)  = 


1  - 


a/c) 


s(fc) 

provided  that  7^  0.  Furthermore, 


(8) 


n(fc+D  =  uW  - 


1  +  O^F  (u^  —  ^4^^  1  i?(fc) 


(9) 


The  scalar  parameter  represents  the  acute  angle  between  the  directions 
generated  by  step  2.1  and  a  chord  step.  Note  that  if  the  directions  s ^ 
and  are  mutually  orthogonal  then  the  chord  step  is  taken  in  full.  The 
composition  of  steps  of  this  kind  has  been  devised  as  a  mechanism  to  generate 


higher-order  nonlinear  algorithms  [17,  24].  In  fact,  if  A ^  =  J ^  and  9 ^  =  1 
for  k  =  0, 1, . . .  then  (9)  becomes 


u(fc+i)  =  u{k) 


-l 


p(k)  _|_  p 


(10) 


for  k  =  0, 1, . . ..  This  recurrence  represents  a  higher-order  modification  of 
Newton’s  method.  Iterates  generated  by  (10)  converge  g-superlinearly  with 
(/-order  3  [24],  These  methods  were  studied  by  Shamanskii  [28]  and  Traub 
[29].  They  pointed  out  that  even  higher-order  methods  can  be  built  out  of 
a  longer  sequence  of  chord  steps  alternated  with  regular  Newton  steps.  In 
a  more  recent  treatment,  Kelley  names  those  methods  after  Shamanskii  and 
compares  the  particular  case  (10)  numerically  against  Newton’s  method  [17]. 

Along  the  lines  of  Gay’s  local  convergence  analysis  for  Broyden’s  method, 
Yang  shows  that  the  NEN  algorithm  converges  n-step  (/-quadratically  for  n- 
dimensional  problems  [34],  Therefore,  as  in  the  linear  case,  the  NEN  method 
converges  twice  as  fast  as  Broyden’s  method.  In  addition,  efficient  limited 
memory  BFGS  implementations  lately  proposed  for  Broyden’s  method  apply 
here  as  well  [22,  30]. 


4  Krylov-Secant  Updates 

Let  A ^  be  an  approximation  of  the  current  Jacobian  J ^  .  We  are  interested 
in  looking  at  a  minimum  change  of  A ^  consistent  with  the  secant  equation 
A (k)sik)  =  pfi+A  _  pik)  and  restricted  to  the  underlying  Krylov  subspace. 
A  basis  for  this  subspace  arises  as  a  result  of  using  a  linear  iterative  solver 
such  as  GMRES  for  solving  the  approximated  Jacobian  system  associated 
with  A™. 

4.1  Secant  Updates  constrained  to  the  Krylov  subspace 

The  Arnoldi  factorization  provides  valuable  information  that  should  not  be 
discarded  at  every  GMRES  restart.  This  observation  has  been  a  cornerstone 
of  several  efforts  to  initiate  the  reuse  or  recycling  of  Krylov  basis  information 
to  generate  improved  restarted  versions  of  GMRES  [5,  23,  26].  We  show 
how  to  reflect  secant  updates  on  the  Jacobian  matrix  without  altering  the 
current  Krylov  basis.  For  the  sake  of  simplicity,  let  us  omit  the  sources 
of  inexactness  induced  by  the  use  of  GMRES  whose  relative  residuals  are 
supposed  to  converge  at  a  predefined  tolerance  (i.e.,  to  a  prescribed  forcing 
term  value). 
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For  a  given  nonlinear  iteration  k,  recall  that  the  Arnoldi  factorization  for 
AW  is  given  by  [27] 

A^VL{k)  =  +  /gjef,  (11) 

where  €  M  holds  the  Krylov  orthogonal  basis  along  its  columns,  Hjk>  £ 
R.  is  an  upper  Hessenberg  matrix  and  /,v  ;  is  a  residual  vector  orthogonal  to 
Vyl\  Equation  (11)  may  also  be  expressed  in  the  more  compact  form: 

d(fc)y/fe)  =  V^Hf\  (12) 


Mk)  = 


K  = 


Consider  the  solution  of  the  following  approximated  Jacobian  equation 
at  the  fc-th  nonlinear  iteration 


A^s{k)  =  -F(k\  (13) 

with  the  GMRES  algorithm.  This  linear  solution  may  be  realized  as  part  of 
an  inexact  Newton  or  Broyden’s  method.  Let  s[k'1  =  s +  vf^y^  be  the 
solution  obtained  after  l  steps  within  the  final  GMRES  restart  cycle.  The 
associated  Krylov  subspace  for  this  problem  is  given  by  K,^  . 

Now,  we  seek  to  use  the  information  gathered  during  the  solution  of  (13)  to 
provide  an  approximation  to  the  system 

Equation  (14)  may  be  seen  as  a  perturbed  version  of  linear  system  (13) 
and  hence,  in  general,  we  can  not  guarantee  that  K^k+1\A^k+1\  rgk+1^)  = 
/cjk\A(k\  r^).  However,  rank-one  updates  to  the  Arnoldi  factorization  of 
(13)  may  be  done  without  destroying  the  Krylov  basis.  That  is, 

(a®  +  vfW  (l-f >)‘)  vf 1  =  if  ’  (flf>  +  zw‘)  +  (15) 

or  equivalently, 

(V^  A^  +  Vl{k)zwt  (V^y  Vt{k)  =  H[k)  +  zw\  (16) 
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for  any  vectors  z,w  G  M1.  Expression  (16)  suggests  a  clearer  way  to  update 
rather  than  A^k\  Note  that  the  current  Jacobian  approximation  ap- 


„(*) 


pears  to  be  updated  by  a  rank-one  matrix  whose  range  lies  on  K,^  ,  r0 

In  order  to  proceed,  it  would  be  convenient  to  express  the  secant  equation 


^(fc+!)s(fc)  _  p(k+ 1)  _  p(k) 


(17) 


in  terms  of  a  solution  strictly  lying  on  the  Krylov  subspace.  In  addition,  if 
the  initial  guess  to  GMRES  is  other  than  zero,  we  may  rewrite  (17)  in  terms 
of  a  system  whose  solution  lies  on  the  Krylov  subspace.  That  is,  equation 
(13)  would  satisfy 


A(fcVfc)  =  - F (fc)  -  A^sik)  =  r{k\ 


(18) 


with  s?''1  =  1 and  thus,  the  corresponding  secant  equation  (17)  would 
become 

A(k+i)s(k)  =  F(k+ 1)  _  rW.  (19) 

Multiplying  both  sides  by  ,  it  follows  that  +  should  satisfy  the 

secant  equation 

(20) 


H, 


(fc+i) 


yW  =  (v^y  +  f3ei, 


where  (3  = 


„(fc) 


Therefore,  the  Krylov  subspace  projected  version  of  the 


secant  equation  (17)  can  be  written  as 

^(vr/*))tF(fc+1)  +  0e i  -  (y^ 


H, 


(k+i) 


=  H[k)  + 


(y  (^') )  ^  y^) 

This  update  in  terms  of  A ^  is  as  follows 


(21) 


^(fc+l)  _  j^{k)  _|_  p 


(a(fc))*s(fc) 


(22) 


with  P  =  Vj<k>  (vfk^  .  We  refer  to  the  update  (21)  as  the  Krylov-secant 
update.  Note  that  the  operator  P ^  is  an  orthogonal  projector  onto  the 
Krylov  subspace  lCm  (^A^,  .  That  is, 


p(k)^  =  p(k)  (Idempotency) . 
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•  ^p(fc)j  =  p(fc)  (Symmetry). 

•  Range^P(fc^  =  Km  (^A^k\  . 


4.2  Krylov  Secant  vs.  Standard  Secant  Update 


It  can  be  shown  that  expression  (22)  corresponds  to  the  solution  of  the 
following  minimization  problem 


min 

BeRnxn 


B  -  A{k) 


F 


subject  to  s(fc) 


p{k)F{k+i)  +rWj 


where  ||-||^,  stands  for  the  Frobenius  norm.  In  other  words,  the  Krylov-secant 
update  produces  the  closest  matrix  in  Frobenius  norm  satisfying  the  secant 
equation  that  lies  on  the  Krylov  subspace  JC^  ^A^k\  r^'j .  This  problem 
has  the  equivalent  formulation  in  terms  of  an  Hessenberg  matrix,  namely, 


min 

GeRmXm 


G  -  H& 


F 


subject  to  Gy  =  F™  +  fa. 


A  geometric  interpretation  may  be  drawn  out  from  this  fact.  Consider 
Q,  the  set  of  matrix  quotients  of  q  =  p(fc+1)  —  F ^  by  s  =  —  s ^ 

defined  by 

Q  =  {B  G  Rnxn  |  Bs  =  y},  (23) 

and  A,  the  set  of  matrices  generating  the  same  Krylov  subspace  /Cm  = 
ICm  (A(fe),rf}).  That  is, 

A  =  {Be  Mnxn  |  ICm(B,  r(0k) )  =  JCm}.  (24) 

The  resulting  matrix  A(k+l'>  in  (22)  may  be  viewed  as  the  nearest  to  A ^  in  X 
among  the  set  of  nearest  matrices  in  X  to  Q.  Furthermore,  if  the  intersection 
of  these  two  sets  is  not  empty,  then  E  X  n  Q.  This  observation  is  key 

to  the  construction  of  least-change  secant  updates  consistent  with  operators 
satisfying  the  standard  secant  condition  and  other  properties  prescribed  by 
a  given  affine  subspace  (e.g.,  sparsity  pattern,  positive  definiteness)  in  Mnx,t 
(see  [7,  8]). 

Figure  2  provides  a  geometric  interpretation  of  Krylov-secant  updates 
and  standard  secant  updates.  Here,  A^-  represents  the  new  approximation 
obtained  by  (22),  Ag  represents  the  new  approximation  obtained  by  (7),  and 
J  is  the  Jacobian  operator  at  the  new  point.  Depending  on  the  behavior  of 
the  function  F.  either  of  these  two  updates  will  yield  an  operator  closer  to 
the  Jacobian. 
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Figure  2:  Geometric  interpretation  of  Krylov-secant  updates  as  compared  to 
standard  secant  updates. 
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4.3  The  Nonlinear  Krylov-Eirola-Nevanlinna  Algorithm 

We  are  now  in  a  position  to  describe  an  inexact  NEN  algorithm  that  ex¬ 
ploits  the  information  left  behind  by  the  GMRES  method.  We  introduce  the 
nonlinear  Krylov-Eirola-Nevanlinna  (nonlinear  KEN)  algorithm  as  follows: 


Algorithm  4.1  (Nonlinear  KEN) 

1.  Give  an  initial  guess  u ^  and  a  Jacobian  approximation 

2.  For  k  =  0, 1, . . .  until  convergence  do 


2.1 


h 


(k) 

ra+l,m’ 


P  = 


2.2 

2.3 


q(k)  = 


H, 


(fc+i) 


(v)(fe))tF(fe+1)  +  /3e1. 

u{k)  |  (<z(fc)-tf,(fc)v<fc))  (i/<*))‘ 
'T :  1  (yW)tj/(fc) 


2.4  Solve 


=GMRES(AW, 


mm 


/Je1  +  i4ifc+1)y 


with  ET,(fc+1)  =  f 


(fc+i) 


Denote  its  solution  by 


2.5  =  V)(fc)y(fc). 

2.6  u(fc+1)  =  +  s(fc). 

2.7  Update  or  construct  Afc+1). 


Note  that  the  basic  difference  between  the  NEN  algorithm  and  the  non¬ 
linear  KEN  algorithm  is  the  way  in  which  the  composed  nonlinear  step  is 
computed.  In  the  nonlinear  KEN  algorithm  there  are  two  updates:  one  in 
n-th  dimensional  space  and  another  restricted  to  the  Krylov  subspace  of  di¬ 
mension  l.  The  computational  implications  are  evident.  Whereas  the  NEN 
algorithm  solves  two,  generally  large,  linear  systems  of  the  same  size,  the 
nonlinear  KEN  solves  one  large  and  one  significantly  smaller  system  (since, 
in  general,  l  <C  n).  For  the  sake  of  updating  the  Jacobian  approximation, 
step  2.7  in  the  nonlinear  KEN  algorithm  may  be  implemented  as  step  2.3  of 
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the  NEN  algorithm.  Another  option  is  to  form  from  scratch  at  the 

new  nonlinear  step  so  that  a  true  Newton  step  is  composed  with  this  Krylov 
secant  update. 

As  one  may  expect  from  the  discussion  of  the  NEN  algorithm  leading 
to  the  update  (9),  it  may  be  shown  that  the  resulting  composite  step  is  the 
Broyden  or  Newton  method  with  a  chord  step  constrained  to  the  Krylov 
basis. 

4.4  A  high-order  Newton— Krylov  algorithm 

An  even  faster  version  of  the  nonlinear  KEN  algorithm  may  be  attained 
by  performing  exhaustive  rank-one  updates  of  the  Hessenberg  matrix  before 
making  the  next  GMRES  call.  The  extent  of  these  updates  is  determined 
by  the  capability  of  the  residual  minimization  problem  to  produce  descent 
directions.  Hence,  we  stop  performing  Hessenberg  updates  when  we  are 
unable  to  generate  a  sufficient  decrease  of  ||.F||.  Details  on  this  procedure 
were  described  in  [33]. 


5  Implementation  issues 

Since  Krylov-secant  updates  imply  several  operations  on  a  restricted  basis 
(or  smaller  space),  computations  may  be  carried  out  in  an  efficient  fash¬ 
ion.  We  proceed  to  mention  the  implications  of  Krylov-secant  updates  for 
preconditioning,  globalization  and  forcing  term  strategies. 

5.1  Preconditioning 

Expression  (16)  shows  that  if  a  preconditioner  is  used,  the  Krylov-secant 
update  (21)  will  reflect  the  implicit  joint  update  of  the  Jacobian  times  its 
preconditioner.  More  precisely,  if  M ^  is  an  operator  whose  inverse  approx¬ 
imates  the  inverse  of  the  Jacobian  (or  an  approximate  Jacobian  operator), 
then  the  Krylov-secant  update  for  right  preconditioning  implies 


A ^  (m™ 


-l 


+  P 


(g(k)_A(fc)s(*))  (#))' 

(«(*))*#) 


(25) 
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with  s ^  =  M^s^k\  On  the  other  hand,  left  preconditioning  reads  as 
follows 


lA^  +  P 


G 


(fc) 


-A^s^)  («(*>)* 

(s(fc))*  g(fc) 


(26) 


Note  that,  in  the  latter  case,  there  are  no  changes  in  the  search  direction 
with  respect  to  the  unpreconditioned  problem.  These  two  implicit  joint 
updates  must  be  treated  with  care  in  order  to  maintain  consistency  in  the 
computations.  The  key  point  is  that,  in  general,  one-rank  updates  do  not 
distribute  with  respect  to  the  product  of  matrices  and  therefore  we  are  unable 
to  guarantee  that  either 


=  A+ (M+)  1  or  +  =  (M+)  1  A+ . 

In  addition,  particularly  for  multiphase  groundwater  flow  purposes,  com¬ 
putation  of  an  updated  version  of  the  preconditioner  may  be  a  time  consum¬ 
ing  task  or  even  impossible  to  carry  out  when  there  is  not  an  explicit  form 
of  such  preconditioner  as  it  happens  in  two-stage  preconditioning  [6,  20] 
multilevel  [4]  or  in  general  inner-outer  iterations  [2],  Assuming  that  the 
preconditioner  is  fixed  after  each  Krylov-secant  update  seems  to  be  the  sim¬ 
plest  and  most  efficient  strategy,  even  though  this  may  introduce  some  extra 
inaccuracies  in  the  computations  [18]. 


5.2  Globalization  Strategy  and  Forcing  Terms 

The  utilization  of  GMRES  offers  the  opportunity  to  avoid  explicit  repre¬ 
sentation  of  the  Jacobian  operator  after  each  Krylov-secant  update.  To 
implement  the  globalization  strategy  we  must  construct  the  interpolating 
polynomial  that  satisfies 

\F'  (0)  =  ^  (F,  Jsm)  =  ( F ,  J  (s0  +  Vmym)) 

=  (F,  -F  -  r0>  +  (F,  Vm+1Hmym )  (27) 

=  -  \\F\\2  -  (F,  Vm+i  (/3e i  -  Hmymy  . 

Hence,  no  explicit  reference  to  the  Jacobian  is  required.  Note  that  (27) 
simplifies  to  ( F,Jsm )  =  —  ||F||2  +  ||rm||2  when  sq  =  0.  Furthermore,  to 
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compute  residuals  associated  with  the  damped  nonlinear  step  we  need  only 
to  compute 


(A(fc)s(fc)) 

2 

p(k) 

2  +  2A(fc) 

+  ( 

A(fc) 

V 

j(k)s(k) 

)2. 

where 


(28) 


2 


(29) 


Finally,  note  that  the  Jacobian  operator  may  also  be  bypassed  in  the 
forcing  term  computation,  since 


p(k) 

- 

(3e i  -  Hmy 

p(k) 

-  (3  <?m+ id 

F 

(*- 1) 

F(k- 

1) 

(30) 


6  Numerical  Experiments 


Now  we  looked  at  a  two-phase  problem  consisting  of  wetting  and  a  non¬ 
wetting  phases.  Primary  variables  are  non-wetting  phase  pressures  and  wa¬ 
ter  concentrations.  The  model  consists  of  a  bottom  hole  pressure  specified 
water  injection  well  at  coordinate  (1,1)  of  the  plane  and  a  bottom  hole  pres¬ 
sure  specified  production  well  in  the  opposite  corner.  All  terms  are  solved 
implicitly. 

Absolute  permeabilities  are  set  according  the  following  diagonal  tensor 
K (in  milidarcies) 


K  = 


100  0  0 

0  500  0 

0  0  50 


Grid  sizes  are  defined  as  DX  =  DY  =  20  ft  and  DZ  =  2  ft.  The  timestep 
increase  factor  is  set  to  1.25  in  order  to  reach  rapidly  the  maximum  DT  size 
of  45  days.  This  allows  the  production  of  appreciable  changes  in  pressures 
and  saturations  and  cause  the  activation  of  the  several  backtracking  steps. 
Table  1  shows  the  rest  of  the  input  parameters  as  specified  in  this  simulation 
study.  Figure  3  show  the  relative  permeability  curves  and  capillary  pressure 
associated  with  the  case  model. 

We  use  the  Gauss-Seidel  two-stage  preconditioner  [18]  for  solving  the 
linear  coupled  system.  In  brief,  this  preconditioner  consists  of  a  decoupling 
stage  defined  by  the  inverse  diagonals  blocks  of  the  coefficients  representing 
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the  contribution  from  oil  pressures,  water  saturation  and  their  interaction 
within  the  same  grid  block.  The  second  stage  consists  of  computing  an 
iterative  solution  of  pressures  and  saturations  from  the  decoupled  system. 
The  solutions  are  combined  in  a  Gauss-Seidel  fashion  and  hence,  the  whole 
procedure  entails  an  inner-outer  iteration  for  the  solution  of  the  original 
coupled  system.  Numerical  experiments  were  carried  out  on  a  Linux-based 
cluster  consisting  of  45  dual  processors  and  interconnected  through  Myrinet 
cards. 


Table  1:  Simulation  input  data. 


Simulation  Time  (days) 

500 

Max.  linear  solver  tolerance 

lxlO-5 

Nonlinear  solver  tolerance 

5.0xl0-3 

Number  of  gridblocks 

48x48x8  =  18, 432 

DT  (days) 

45 

Water  compressibility  (psi-1) 

3.3xl0-6 

Non- wetting  compressibility  (psi-1) 

4.2xl0-5 

Water  density  (lb/ft3) 

62.00 

Non- wetting  density  (lb/ft3) 

48.00 

Water  viscosity  (cp) 

0.23 

Non-wetting  viscosity  (cp) 

1.6 

Initial  water  saturation 

.2 

Initial  pressure  (psi) 

1800 

Injector  BHP  (psi) 

2000 

Producer  BHP  (psi) 

1600 

Table  2  shows  how  in  this  relatively  simple  problem,  the  Krylov-secant  is 
able  to  produce  important  reductions  in  the  computation.  This  is  reflected 
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Table  2:  Performance  comparison  between  Newton  method  and  the  nonlinear 
KEN  algorithm. 


Method 

Linear  Iter. 

Nonlinear  Iter. 

K-S  updates 

Back’s 

CPU  time  (s) 

Newton 

160 

1076 

- 

2 

779 

N.  KEN 

107 

780 

92 

6 

482 

in  the  nonlinear  KEN  algorithm  as  a  smaller  number  of  nonlinear  itera¬ 
tions,  linear  iterations  and  CPU  time  in  constrast  to  the  Newton  method.  It 
is  worth  mentioning  that  the  Krylov-secant  steps  produce  marginally  more 
backtracking  steps  than  Newton  for  the  successful  iterations.  However,  New¬ 
ton  produces  several  backtracking  steps  at  the  beginning  of  the  simulation 
that  causes  a  halving  of  the  timestep.  There  were  indications  that  the  non¬ 
linear  KEN  may  have  a  larger  convergence  region  than  Newton  under  certain 
situations.  This  fact  has  motivated  other  researchers  to  combine  Broyden 
or  Picard  iterations  with  the  Newton  method  in  order  to  develop  efficient 
continuation  methods  for  timestepping. 

Figures  4  compares  the  performance  of  both  methods  throughout  the 
whole  simulation.  Clearly,  the  nonlinear  KEN  outperforms  Newton’s  method. 
We  can  observe  how  the  difference  in  performance  builds  up  as  the  simula¬ 
tion  progresses.  As  one  might  expect,  the  CPU  time  response  is  a  direct 
consequence  of  the  number  of  linear  and  nonlinear  iterations. 

7  Conclusions  and  Further  Work 

We  have  proposed  a  new  family  of  Krylov-secant  updates  with  the  potential 
for  a  high-order  of  convergence  in  environmental  problems.  These  updates 
do  not  explicitly  rely  on  the  computation  or  use  of  the  Jacobian  matrix  that, 
in  general,  involves  millions  of  floating  point  operations.  We  described  the 
nonlinear  KEN  algorithm  and  a  higher  order  version  of  it  that  showed  to  be 
superior  to  Newton’s  method  for  the  numerical  cases  explored  in  this  paper. 
The  nonlinear  KEN  seems  to  be  worthy  of  consideration  in  fully  implicit 
formulations  where  the  solution  of  a  nonlinear  algebraic  system  must  be 
efficiently  overcome. 

Nevertheless,  the  potential  for  improving  the  present  method  is  enor¬ 
mous.  More  robust  versions  may  be  developed  by  enriching  or  refreshing  the 
Krylov  basis  in  each  of  the  nonlinear  iterations.  Advances  along  these  lines 
have  been  used  to  iteratively  solve  linear  systems  with  different  right  hand 
sides  or  a  sequence  of  slowly  changing  linear  systems  in  the  matrix  and  the 
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Performance  in  Terms  of  Linear  Iterations 


Performance  in  Terms  of  Nonlinear  Iterations 


Performance  in  Terms  of  CPU  time 


Figure  4:  Comparison  of  linear  (top),  nonlinear  iterations  (middle)  and  CPU 
time  (bottom)  between  the  nonlinear  KEN  and  Newton  method. 
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right  hand  side  [5,  23,  26].  We  also  find  promising  the  use  of  deflation  strate¬ 
gies  oriented  to  capture  the  physics  of  the  problem  through  eigen-information 
[1,  12]. 

We  remark  that  our  results  are  promising  but  must  be  evaluated  un¬ 
der  more  stringent  physical  situations  and  at  a  larger  scale.  Among  tar¬ 
get  applications  we  consider  useful  the  insights  coming  from  two-phase  air- 
water,  chemical  processes  and  their  coupling  to  geomechanics.  The  different 
strengths  in  nonlinearities  between  the  typical  coupled  equations  of  the  sub¬ 
surface  model  (i.e.,  one  for  the  pressure  potential  or  flow  driving  force  and 
another  or  others  for  the  saturations/  concentrations)  and  the  coupling  itself 
between  reservoir  variables  are  sufficiently  good  reasons  to  investigate  more 
of  the  preconditioning  ideas  discussed  in  this  project. 
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